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Abstract 

In a variety of neutron lifetime experiments, in addition to f3— decay, neutrons 
can be lost by other mechanisms including wall losses. Failure to account for 
these other loss mechanisms produces systematic measurement error and as¬ 
sociated systematic uncertainties in neutron lifetime measurements. In this 
work, we construct a competing risks survival analysis model to account for 
losses due to the joint effect of /3 —decay losses, wall losses of marginally 
trapped neutrons, and an additional absorption mechanism. We determine 
the survival probability function associated with the wall loss mechanism by 
a Monte Carlo method. We track marginally trapped neutrons with a sym- 
plectic integration method that assumes neutrons are classical particles. We 
model wall loss probabilities of Ultracold Neutrons (UCNs) that collide with 
trap boundaries with a quantum optical model. Based on a fit of the com¬ 
peting risks model to a subset of the NIST experimental data, we determine 
the mean lifetime of trapped neutrons to be approximately 700 s - consid- 
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erably less than the current best estimate of (880.1 ± 1.1) s promulgated 
by the Particle Data Group [1]. Currently, experimental studies are under¬ 
way to determine if this discrepancy can be explained by neutron capture 
by 3 He impurities in the trapping volume. We also quantify uncertainties 
associated with Monte Carlo sampling variability and imperfect knowledge 
of physical models for neutron interactions with materials at the walls of 
trap as well as beam divergence effects. Finally, in a Monte Carlo experi¬ 
ment, we demonstrate that when the trapping potential is ramped down and 
then back up again, systematic error due to wall losses of marginally trapped 
neutrons can be suppressed to a very low level. Monte Carlo simulation stud¬ 
ies indicate that this ramping strategy is more efficient than an alternative 
ramping scheme where UCNs are produced when the held is fully ramped, 
and then increased to its maximum value after the trap is filled. The survival 
probability formalism developed here should be applicable to other experi¬ 
ments where neutron (or other particle) loss mechanisms are non-trivial (i.e., 
where the associated survival probability function with the loss mechanism 
is non-exponential). 

Keywords: Chaotic scattering 
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Symplectic integration 

PACS: code 81.07.Gf, 07.05.Kf, 07.79.-v, 07.79.Lh 


1. Introduction 

According to the Particle Data Group, the current best estimate of the 
mean lifetime of the neutron r n is (880.1 ± 1.1) s [1]. This uncertainty is 
largely due to systematic effects. Hence, reduction of systematic uncertain¬ 
ties is key to reducing the overall uncertainty of the mean lifetime of the 
neutron [2], To emphasize this point, we note that the difference between 
neutron lifetime estimates determined by beam and bottle experiments is 
(8.4 ± 2.2) s. This difference is likely due to unresolved systematic errors [3]. 
In trapping experiments, marginally trapped neutrons can escape a trapping 
volume before they fd— decay. Failure to account for this loss mechanism 
leads to systematic measurement error and associated systematic uncertain¬ 
ties in neutron lifetime measurements [4-9]. In this work, we focus on wall 
losses of marginally trapped Ultracold Neutrons (UCNs) in a magnetic trap- 
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ping experiment at NIST. Since the probability of losing a neutron when 
it interacts with materials at a wall depends, in part, on neutron energy, 
the survival probability of a marginally trapped UCN is non-exponential. 
We estimate the systematic error in the neutron lifetime estimate due to 
marginally trapped UCNs for a particular subset of the NIST data where 
the trapping potential is static. We stress that systematic error (and as¬ 
sociated systematic uncertainty) associated with marginally trapped UNCs 
can be reduced by field ramping strategies to very low levels with the cost 
of increasing random uncertainty. Hence, the systematic error we report for 
the static potential case does not apply to cases where marginally trapped 
UCNs are purged by field ramping methods. 

In the NIST experiment [5,6] an UCN is produced during a trap-filling 
stage when a 12 K neutron is scattered to near rest in liquid helium by 
single phonon emission [10]. After filling the trap and blocking the cold 
beam that produces UCNs, light produced by neutron decay and background 
events is detected by a pair of phototubes. A marginally trapped (or above 
threshold) neutron is an UCN with sufficient energy to escape the trap by 
interacting with materials at the boundary of the trap [11,12,13]. A UCN 
with insufficient energy to escape the trap by the wall loss mechanism is 
a below threshold UCN. Since a marginally trapped UCN can be lost to 
the walls before it f3— decays, the decay rate of marginally trapped UCNs is 
not exponential. Hence, when wall losses are non-negligible, on average, the 
observed /3 —decay rate of neutrons is also non-exponential. Thus, if one were 
to fit an exponential model to such decay event data, one would expect the 
associated estimate of the neutron lifetime to be biased. 

Here, we develop a competing risks survival analysis model [14,15,16,17] 
that accounts for loss of UCNs due to (3— decay, wall losses, upscattering 
and other absorption mechanisms. Determination of the survival probabil¬ 
ity function associated with the wall loss mechanism is the key problem 
addressed by this work. We determine the survival probability function as¬ 
sociated with wall losses by a Monte Carlo method as described in [18]. 

Based on the conditional survival probability function of an UCN that 
survives the filling stage, we construct a prediction model for observed data. 
If /3 —decay and wall losses were the only loss mechanisms, and we had per¬ 
fect knowledge of the survival probability of marginally trapped UCNs as 
well as possible backgrounds, a fit of this model to data would yield a nearly 
unbiased estimate of the true neutron lifetime given a sufficient amount of 
data. In practice, there are other loss mechanisms including upscattering 
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and neutron absorption by 3 He that affect the mean lifetime of neutrons in 
the trap. Still, given perfect knowledge of the survival probability functions 
associated with these loss mechanisms, one could correct the measured neu¬ 
tron lifetime of trapped UCNs and recover the true neutron lifetime in an 
ideal noise-free experiment. We analyze a subset of the NIST experimental 
data where neutrons are confined in a static potential to quantify various 
sources of systematic measurement error and associated uncertainty. For the 
data analyzed, the measured lifetime of trapped UCNs is approximately 700 
s. We attribute the discrepancy between our estimate and the current best 
estimate of 880.1 s to an additional absorption loss mechanism. Currently, 
experiments are underway to determine if 3 He impurities are sufficiently high 
to explain the discrepancy. 

In some runs of the NIST experiment, we vary the trapping potential to 
purge marginally trapped UCNs. In a simulation experiment, we demon¬ 
strate that held ramping can in principle reduce systematic measurement 
error and associated uncertainty due to marginally trapped UCNs to very 
low levels at the expense of purging below threshold UCNs. For a particular 
example, we demonstrate that the above ramping strategy is more efficient 
than one where neutrons are trapped when the held is fully ramped and 
then increased to its maximum value. By more efficient, we mean that the 
expected number of below threshold UCNs after ramping is higher for our 
ramping strategy as compared to the alternative strategy given that both 
methods purge marginally trapped UCNs with nearly the same efficiency. 

In this paper, we hrst present our physical and computational methods 
for simulating UCN trajectories and loss probabilities of UCNs due to wall 
interactions. Next, we develop our competing risks survival analysis model 
and apply it to experimental data acquired at NIST. We also quantify un¬ 
certainties associated with Monte Carlo sampling variability and imperfect 
knowledge of physical models for neutron interactions with materials at the 
walls of trap and beam divergence effects. Finally, in a Monte Carlo simu¬ 
lation experiment, we demonstrate that held ramping can reduce systematic 
error due to wall losses of UCN to a very low level at the expense of reducing 
the expected number of trapped UCNs after ramping. 


4 



2. Physical model 

2.1. Trapping potential model 

In the NIST experiment, an UCN with sufficiently low energy is trapped 
in the potential field produced by the interaction of the magnetic moment of 
a neutron and a spatially varying magnetic field (Figure 1) and gravity. For 
tracking UCNs in the trapping volume, assuming adiabatic spin transport, 
the potential is 

V(x) =//|B(x)| + m n gy, (1) 

where B is the magnetic field, m n and /i are the mass and magnetic moment 
of the neutron respectively, x = (x, y, z ) is spatial location of the neutron in 
the Cartesian coordinates shown in Figure 1, and g is the acceleration due 
to gravity. 

We define a nominal trapping volume corresponding to — z Q < z < z 0 
where z is the axial coordinate and z 0 = 37.5 cm (see Figure 1). An UCN 
with total energy (kinetic plus potential) greater than the minimum of the 
potential Win on the boundary of the nominal trapping volume is defined to 
be a marginally trapped (or above threshold) UCN. An UCN with energy 
less than or equal to Win is a below threshold UCN. We stress that there 
is not a physical boundary at z = -37.5 cm. However, as shown in Figure 
1, as z is reduced from z ~ -37.5 cm to -60 cm, the trapping potential 
decreases dramatically. For modeling purposes and to speed up Monte Carlo 
simulations, it is important to define a nominal trapping volume so that all 
UCN with initial energy less than the minimum value trapping potential on 
the boundary of the nominal trapping volume never interact with materials 
at the ZZ at z = -60 cm and z = 37.5 cm, nor the cylindrical boundary r = 
5.6 cm. If we choose the lower axial boundary to be much less than -37.5 cm, 
this condition is violated. Different choices of the minimum axial location 
of the nominal trapping volume affect results. In Section 4.3, we quantify 
this effect. For a typical static trap configuration of the NIST experiment, 
Win = 139 neV. We assume that probability distribution function of the 
initial speed |u| of an UCN produced in the trap has a quadratic form [18,19] 
f(\v\) oc |u| 2 for low velocities of interest. 

We determine a neutron trajectory based on its initial position and mo¬ 
mentum, by solving the classical equations of motion 

P = ^(x) = —VV(x), (2) 
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and 


P 


( 3 ) 



with an optimal fourth order symplectic integration scheme [18,19,20,21], 
We predict |B| at arbitrary points in the trapping volume with a three- 
dimensional tensor-product spline interpolant [22] where the order of the 
spline is four in each direction. We determine the tensor-product B-spline 
coefficients from values of |B| computed on a grid by a code that solves the 
Biot-Savart law numerically based on the geometry of the solenoid and cur¬ 
rent bars that produce the magnetic held. With this tensor-product method, 
we evaluate the potential and its gradient at arbitrary locations within the 
trap. 

In some experiments, we ramp the quadrupole held B^ while keeping the 
solenoid held B s constant. Given that the filling stage ends at time the 
ramping factor for the quadrupole held is R(t — U.), the magnitude of the 
B-held varies as 


|B (t -t L )| = |B S + R(t - t L )B„\, (4) 

where B q>max is the maximum value of the quadrupole held. For ramping 
cases, we develop a tensor-spline method to model the gradient of |B(t — £l)| 
based on tensor-spline models for each component of B A and B, hmax . 

2.2. Wall loss model 

When an above threshold UCN collides with the cylindrical wall or endcap 
boundaries, we assign it a loss probability pi oss according to a model based 
on assumed material properties [11,12,13]. After k collisions, the empirical 
survival probability of the UCN is 

k 

PsVLIv(k) = jjl Ploss(^)] - (3) 

i= 1 

We track each above threshold UCN until its empirical survival probability 
drops below 10” 9 . Due to chaotic scattering effects, the symplectic integra¬ 
tion prediction for a trajectory of an UCN does not converge in general as 
the time step parameter in the integration code is reduced [19]. Here, we 
assume that the mean survival probability at any time t for an ensemble of 
UCNs does not depend on the time step parameter even though the pre¬ 
dicted survival probability of a particular UCN may depend on the time step 
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parameter. Although we are unaware of any proof that this assumption is 
true, it seems reasonable. For more discussion of this point, see [19]. 

We model the interaction of the neutron with the materials on the trap 
boundaries with a one-dimensional optical model based on Schrodinger’s 
equation. In this approach, we assume that each neutron energy is suffi¬ 
ciently low so that its wavelength is very large compared to the spacing 
between nuclei in the wall materials. Hence, coherent effects are significant 
and interactions are well predicted by an optical model where the neutron 
potential is V — iW. For the cylindrical walls, we model the the neutron 
potential as due to layers of different homogenous materials following [11], 
The materials at the endcaps at z = —60 cm and z = 37.5 cm are Teflon 
FEP Q and acrylic respectively. Given that a UCN with velocity v crosses 
the trap boundary at location x and that the surface normal for the trap 
boundary at x is n, we define E± = \m n |v ■ n| J . For the endcaps, we model 
the probability of diffuse reflection off the wall based on Eq. 2.71 of [11] as 

_ E_ L - y/EjJ 2q - 2(K - E _ L )) + a 
Pscat ~ E ± + ^/E ± (2a - 2(14 - E ± )) + a 


where 


«= v/(K - E ± y + w 1 2 , 

and K = V — Vff ( , where Vn e = 15.98 neV. For the Teflon material, V = 27.8 
neV and W = 1.39e-04 neV. For the acrylic material, V = 121.04 neV and 
W = 4.74e-05 neV. 

The material that coats the cylindrical walls of the trap is modeled as mul¬ 
tilayer of tetraphenyl butadiene (TPB) (C 28 H 22 ), Gore-Tex (C 2 F 4 ), graphite 
(C), and boron nitride (BN). The TPB used in the experiment is not deuter- 
ated and therefore contains a substantial amount of hydrogen, which lias a 
very large incoherent scattering cross section. Since the optical model does 
not account for incoherent scattering, we add an additional term to account 
for it. The real and imaginary components of the augmented potential for 


1 Certain materials are identified in this paper to foster understanding. Such iden¬ 
tification does not imply recommendation or endorsement by the National Institute of 

Standards and Technology, nor does it imply that the materials identified are necessarily 
the best available for the purpose. 
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where for the ith nuclear isotope, N t is nucleus density; a* is the coherent 
scattering length; u* bs is the absorption cross section; and cr} oss accounts for 
incoherent scattering losses. Since approximately half of the neutrons that 
undergo incoherent scattering will be scattered back into the 4 He and the rest 
will be lost, c i oss is set to half of the estimated total incoherent cross section. 
After this modification, we solve the appropriate differential equation and 
determine the loss probability of the neutron as a function of Ej_ (Figure 
2). Our model does not include the possibility of surface contamination [23] 
as we do not currently have a way of characterizing surface contamination 
quantitatively. However, for the most part, such contamination would lead 
to an additional marginally-trapped neutron loss mechanism. We expect 
that such an additional loss mechanism would reduce the systematic effect 
associated with marginally trapped neutrons. 

After each wall collision, we scatter the neutron back into the trapping 
volume. Since the surface of the walls that scatter neutrons is rough, we 
model the reflection of the neutron with a Lambertian model [24] rather 
than a specular model that is appropriate for a perfectly smooth surface. 
In optical applications of the Lambertian model, the intensity of a reflected 
signal is proportional to cos(0) where 9 is the inclination angle between the 
surface normal of the emitter and the direction of the reflected radiation. 
Hence, in our simulation studies, it follows that the cumulative distribution 
function (CDF) for the inclination angle between the surface normal of the 
wall and the direction of reflected neutron 9 is 


Fl(6) 


1 — cos(2 9) 
2 


where 0 < 9 < |. The azimuthal angle (j) is a uniformly distributed random 
variable between 0 and 2n. In contrast, for a diffuse reflection model based 
on a uniform intensity model, the CDF of the direction cosine of the reflected 
neutron is 


Fir(cos(9 )) = cos($), 



where 0 < cos(d) < 1. Later, in Section 4.3, we quantify the variation of 
the estimated lifetime of trapped neutrons for different neutron reflectivity 
models. 


3. Survival Analysis Model 

3.1. Mathematical preliminaries 

In survival analysis, for any loss mechanism, the loss time T of an object 
(in our case a neutron) created at time t = 0 is a random variable with 
survival probability S(t) where 

S(t) = Pr (T > t). 


For a continuous S(t), one can define the hazard function A (t) which is the 
instantaneous loss rate at time t given that the object of interest survives 
until time t as follows. 


m 


lim 

Ai-s>0 


Pr(t < T <t + At\T > t ) 
A t 


where Pr (t < T < t + At\T > t) is the conditional probability that the loss 
time T falls in the interval [t, t + At] given that T is no less than t. Based on 
the well known conditional probability equality Pr(A|U) = Yv(AC\B) / Pr(5), 
one gets the following well known expression for the hazard function 

m = * iim + At > = -dhsM. 

S(t) At->o At d t 

For the neutron / 3 — decay loss mechanism, the associated hazard function 
is \/3 = A- where r n is the neutron lifetime. For the mechanism associated 
with wall losses of marginally trapped UCNs, we expect X(t) to vary with 
time since high energy UCNs should, on average, be lost at a higher rate 
compared to low energy UCNs. In general, for any loss mechanism, 


S(t) = exp(—A (t)), 
where the cumulative hazard function A(f) is 

A (t) — f X(t)dt. 

Jt =o 
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Thus, the survival probability function associated with /3 —decay is Sp(t) = 
exp(-A /3 r n ) = exp(-^r-). 

In a competing risks model, the multivariate survival probability function 
is S(ti, t 2 , ■ ■ ■, tx) where 


S(t\, t 2 , • • •, tx) — Pr(Ti >ti,T 2 > t 2 , ■ ■ ■ Tx > tx)- 


The term T t is a random variable representing the loss time associated with 
the itli loss mechanism. In general, for i ^ j, the random variables T % and 
Tj need not be independent. The actual loss time of the object of interest is 
T = min(Ti, T 2 , ■ ■ ■ Tx) where K is the number of loss mechanisms. One can 
recover the survival probability function for the jtli loss mechanism at time 
t by evaluating the multivariate one at tj = t and ti = 0 for all other i ^ j. 
The cause-specific hazard function for the jth loss mechanism is A j where 


A jit) 


d log S 


dtj 


tl=t2=—,tK=t 


(7) 


For the NIST experiment, we assume that all neutron loss mechanisms are 
independent. Given this independence assumption, we can write S(t\, t 2 , * • •, tx) 
Y\U S,j{t,j) where Sj(t) is the survival probability function associated with 
the jth loss mechanism. Further, for each independent loss mechanism, the 
cause-specific loss mechanism for the jth loss mechanism is 


A jit) 


dlogSj(t) 

dt 


( 8 ) 


For our problem, the joint hazard and survival probability functions of an 
UCN due to all loss mechanisms at time t are A j(t) and exp(— , A j(t)). 

respectively. 


3.2. Conditional survival probablity 

If an UCN is created at time t — 0, the conditional survival probability 
associated with the first loss mechanism given the the UCN survives all loss 
mechanisms until at least time t l is 


S(t,t L ,t L , ■ ■ ■ t L ) / S(t L ,t L ,t L , ■ ■ -t L ), 

for t>ti,. In the NIST experiment, we model the creation time of any UCN 
during the filling stage as a uniform random variable between t — 0 and 
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t — t^ where tj j is the time spent hlling (loading) the trap. For the hrst loss 
mechanism, the conditional survival probability for an UCN that survives 
the hlling stage is 

Si(t\T >t L ) = [j S(t — s,t L — — s)ds] / 

Js= o 


[f S(t L - s,t L - s, • • • ,t L - s)ds] (9) 

Js =0 

A similar statement applies to the other loss mechanisms. Aside from [3— decay 
and wall losses, neutrons can be lost by absorption processes associated with 
impurities (primarily 3 He capture ) and by upscattering. We expect the haz¬ 
ard function associated with upscattering to vary with temperature. How¬ 
ever, for the ideal case where temperature is constant, the hazard function 
for upscattering is a constant \ u = —. From hrst principle arguments, the 

'Tu 

hazard function associated with an energy-independent absorption process is 
a constant A a = — . We define A* to be sum of the three hazard functions 

3~a 

associated with /3— decay, upscattering and an additional energy-independent 
absorption process. Each of these three hazard functions is a constant in our 
model. Thus, 


A* 


1 

T* 


1 

— + 
r n 


1 1 

-1-• 

A A 


( 10 ) 


We stress that A* does not account for the hazard function associated with 
the wall loss mechanism. Next, we develop a model that enables us to directly 
estimate A* from experimental data given that we have an estimate for the 
survival probability function S'm(^) associated with the wall loss mechanism. 

Following arguments in [18], given the production rate of below and above 
threshold UCNs during the hlling stage are r_ and r + , we predict the ex¬ 
pected number of below threshold UCNs at the end of hlling stage (f^) as 

(■ N_(t L )) = r_f exp(A*(s — t L ))ds. (11) 

Js =0 

The predicted number of above threshold UCNs at the end of the hlling stage 
is 
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(N + (t L )) = r ' + f S M (t L - s) exp(A*(s - t L ))ds, (12) 

Js =0 

where Sm( t) is survival probability function associated with wall losses of 
above threshold UCNs. As stated earlier, S'm(^) has a non-exponential form 
in general. 

Given that an above threshold UCN survives the filling stage, the condi¬ 
tional survival probability associated with wall losses is 

r tL 

S^(t) = S M (t\T >t L ) = [ S M (t ~ s ) exp(A*(s - t L ))ds\ / 

Js =0 

[f S M (t L - s) exp(A*(s - t L ))ds]. (13) 

.7s=o 

For loss mechanisms with exponential survival probability functions, S(t\T > 
tL ) = S(t — ii). Hence, for times t > tL , 

m)) = (N.(t L )) [ 1 + 5 +( t ) | exp(-A.(« - ij). (14) 

We can rewrite the above as 

(N(t)) = {N_(t L )) [ 1 + A (t)] exp(-A,*(t - t L )), (15) 

where the time-dependent “distortion 1 ' term A is 

A(t > = s " (t) ' (16) 

Based on the above, the expected loss rate of neutrons lost due to /3— decay, 
absorption, and npscattering are rp(t ) = ,r a (t) = , and r u (t) = 

respectively. Given that upscattering losses are unobservable and ab¬ 
sorption events yield events that are indistinguishable from (3 —decay events, 
the overall predicted detection rate is 

r*.(*) = (N(t)) P- + (17) 

where pp and p a are detection probabilities for two loss mechanisms. 
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3.3. Contamination ratio 

Following [18], we can decompose (. N(t )) into exponential and non-exponential 
components as follows. 


(N(t)) = fexpit) + f c (t), (18) 

where 

f expit) = {N_(t L )) [1 + Attend)] exp(-A.* (t - t L )), (19) 

and 

feit) = (N_(t L )) [A(f) - A(f end )] exp(—A*(f - t L )). (20) 

The ratio of the non-exponential and exponential terms can be regarded as 
a contamination ratio r c where 

A{t) — A (t en d) 
c[) 1 + A (t end ) 

We emphasize that A (t) and f c (t) are nonlinear functions of A*. In Sec¬ 
tion 4, we estimate A* directly from experimental data given a Monte Carlo 
estimate of S'm(^) based on a physical model for the wall loss probability and 
how surviving neutrons are reflected back into the trapping volume. The 
uncertainty of A* depends, in part, on imperfect knowledge of: the wall loss 
probability model; how surviving neutron reflect off the walls; the spatially 
varying fluence of the thermal beam that produces UCN in the trap; and the 
usual counting statistics variability in the observed data. 

4. Experimental Application 

4-1. Prediction model 

For experimental data corresponding to a particular subset of runs from 
the NIST experiment, we estimate A* by fitting a model to experimental data 
based on Eqns. 15 and 17. The primary (background plus neutron events) 
data are acquired for a static trapping potential. Background measurements 
are also acquired in non-trapping runs and subtracted from the primary ob¬ 
servations [5]. Our primary goal is to understand systematic measurement 
errors and associated uncertainties associated with marginally trapped neu¬ 
trons and other loss mechanisms. Hence, we analyze data corresponding to 
an experiment where we did not ramp the magnetic held in order to maximize 
the systematic effects of marginally trapped neutrons. 


( 21 ) 
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More specifically, for bins that are 1 s wide, we predict the number of 
background-corrected counts in the kth bin as 


n k = AXJ t [ 1 + A (t k )] exp(-A*(4 - t L )) + r bg 5 t , (22) 


where St — 1 s is the resolution at which we determine survival probabilities 
by our Monte Carlo method, t k is the midpoint of the kth time bin, and 
A, r bg and A* are adjustable model parameters determined by our modeling 
fitting procedure. Since the width of the bins for the observed data is 15 s, 
we aggregate predictions at the 1 s scale to get predictions at the 15 s scale 
of interest. 

Given that we estimate these parameters to be A and A, we can predict 
the expected number of below threshold trapped UCNs at t — t k as 


(N-(t L )) 


AX 


tot 


Pp(\ot — S a — X u ) + Pa X a , 


(23) 


where the hazard functions A a and X u and detection probabilities p a and pp 
are determined from other experiments and/or theoretical arguments. 


4-2. Estimation details 

In the NIST experiment, the cold neutron beam that produces the UCNs 
is collimated. Based on an uncollimated measured beam profile and knowl¬ 
edge of the geometry of the collimator, we estimate a spatially varying neu¬ 
tron fluence image at the entrance to the detector (Figure 3) and an asso¬ 
ciated probability density function for the intersection of any neutron tra¬ 
jectory and the plane at z = —60 cm. Based on this probability density 
function, we simulate intersection points with the Von Neuman rejection 
sampling method [25,26]. For each intersection point, we simulate a neutron 
velocity direction vector v that has length 1. For the case where there is no 
beam divergence, v is parallel to the axial direction of the trap. For the case 
of non-zero divergence, we simulate v such that its direction cosine with re¬ 
spect to the ^-direction is uniformly distributed in the interval (cos(# max ), 1). 
That is, simulated realization of v fall within a cone. The location of the 
UCN produced by a neutron is (x 0 ,y 0 ,z 0 ) + L sim v where (x 0 ,y 0 ,z 0 ) is the 
simulated location of the neutron at z 0 = -60 cm, and L sim is the simulated 
distance traveled by the neutron before a UCN is produced. Given the ini¬ 
tial location of the UCN, we simulate its initial velocity as described earlier. 
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Given the initial velocity and position of an UCN, we determine its trajectory 
with the symplectic integration method described in Section 2.1. 

For the case of no beam divergence, the production rate of above threshold 
UCNs, r + (Eq. 12), produced in the trap at random times in the interval 
(0 ,t L ) is 4.2 times larger than the production rate, r_ (Eq. 11), of below 
threshold UCNs for the nominal trapping volume defined by -37.5 cm < 
37.5 cm (Figure 4). The energy range of simulated above threshold UCNs is 
139 neV to 246 neV. In general, the UCN wall collision rate increases with 
energy (Figure 5). Based on Monte Carlo estimates of S'm(^) at discrete 
times (0,1,2, • • •, 5500) s, we determine A (t) (Eq. 16) on a discrete time 
grid (Figure 6). Recall, the theoretical value of A {tif) equals • Even 

though the relative production of above and below threshold UCNs is 4.2, 
(7V + (qjj approximately 0.41. That is, a large fraction of above threshold 
UCN is lost during the filling stage. Since UCNs with energy greater than 
246 neV would be lost to the walls in approximately a few seconds or less, 
such high energy UCNs would have a negligible effect on A (t) for £ — > 10 

s. Thus, extending the maximum energy of simulated above threshold UCN 
would have a negligible effect on our estimate of A*. 

We average 40 independent estimates of Smit) from independent Monte 
Carlo experiments to get an overall estimate of S'^'(t) and A(f). Given a wall 
loss probability model that accounts for incoherent scattering, we determine 
the mean lifetime of the neutron in the trap to be 700 s with a standard 
uncertainty of 57 s (Table 1, Figure 7). In this standard approach, the model 
is assumed to be valid and deviations between observations and predicted 
values based on perfect knowledge of the model parameters are clue to count¬ 
ing statistics. To quantify the component of uncertainty due to imperfect 
knowledge of .SA/(£) due to sampling variability, we simulated bootstrap [27] 
replications of our estimate of 5 'm(^) by resampling with replacement the 
40 independent estimates of S'm(^) that we averaged in the first study. For 
each bootstrap replication of the average 5 'm(^), we refit our model to the 
same observed data. The bootstrap estimate of uncertainty due to sampling 
variability of Sm (£) is 2.4 s. 

4-3. Systematic effects 

For comparison, when we set A(t) = 0, i.e., when we neglect wall losses, 
we estimate the lifetime to be 655 s with an estimated uncertainty of 51 s 
(Table 1). When we estimate S'm(^) based on a wall loss model that neglects 
incoherent scattering effect for the cylindrical boundary, we estimate the 
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neutron lifetime to be 708 s with an associated 1-sigma uncertainty of 59 
s. Based on this, we estimate the component of uncertainty contributed by 
imperfect knowledge of the wall loss probability model to be 8 s. 

For the primary wall loss model that accounts for incoherent scattering, 
we estimate A* for different definitions of the trapping volume. In particular, 
we vary the lower axial boundary of the trapping volume, z mmj from -37.5 
cm to -35 cm and -40 cm (Table 2). We fit a linear model to predict the 
expected value of A as a function of z nnn . The estimated slope of this linear 
model and its associated uncertainty are -0.55 s/cm and 0.40 s/cm. 

In a similar study, we vary the time-step parameter in the symplectic 
integration algorithm (Table 2). The estimated slope of a linear model to 
predict the lifetime as a function of the time step is 2.9 xlCT 4 where the 
uncertainty is 3.79 xlO -4 . Based on this result, the expected difference in a 
lifetime estimate based on a simulation where the time step parameter d t = 
5.0xlO -5 and where df = e where e is arbitrarily close to 0 s, is 1.4 (1.9) s. 

Variation in the beam divergence parameter 9 max produces a statistically 
significant variation in the lifetime estimate (Figure 8) because the distribu¬ 
tion of the initial potential energy of an UCN created in the trap depends 
on 9 max . Thus, even though the distribution of the initial kinetic energy of 
an UCN produced in the trap does not depend on 9 max , the distribution of 
the initial total energy of the UCN depends on 9 max . We estimate 9 max to 
be 3.0 degrees by matching the expected value of 9 in the simulation to that 
predicted based on analysis of the dispersion of scattering plane orientations 
in mosaic crystals relative to their mean value [28]. In [28], the angle between 
the random orientation of a particular scattering plane and the mean orien¬ 
tation is a truncated Gaussian. Since the slope of the fitted line in Figure 8 is 
—0.97(0.16) s deg -1 , we estimate the systematic error due to ignoring beam 
divergence effects to be approximately 2.9 s. As a caveat, to estimate an 
uncertainty, we assume a one-to-one relationship between the standard devi¬ 
ation of the random mosaic crystal orientations and the standard deviation 
of direction cosines in our simulation study. Further, in our beam divergence 
study, the probability distribution function for the direction cosine of a UCN 
at z = -60 cm is a uniform distribution. In contrast, for the mosaic crystal 
model, the angle between the mean orientation and any random orientation 
is a truncated Gaussian distribution. 

Our estimate of A* systematically depends on the assumed model for how 
surviving neutrons are reflected back into the trapping volume (Table 3). 
Because of surface roughness effects, the two diffuse models are much more 
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Table 1. Model parameter estimates for NIST experimental data for a static 
potential. Survival probability function associated with wall loss determined 
for two scenarios; incoherent scattering in neutrons in materials is either 
neglected or accounted for. Beam divergence neglected. The time step pa¬ 
rameter for all cases is dt = 5.0e-05 s. 


case 

account for 

account for 

A 

T * (s) 

n g (s x ) 

X 2 /df 

p-value 


incoherent 

marginally 






scattering 

trapped neutrons 






a 

yes 

yes 

57943 (3883) 

700 (57) 

1.8 (1.8) 

1.080 

0.2234 

b 

yes 

no 

76114 (4531) 

655 (51) 

2.3 (1.7) 

1.085 

0.2118 

c 

no 

yes 

41558 (2856) 

708 (59) 

1.8 (1.8) 

1.082 

0.2187 


plausible than than the specular model. Hence, the difference between the 
estimates for the diffuse models ( 1.3 s with an uncertainty of 1.5 s ) is relevant 
for estimation of a systematic uncertainty due to imperfect knowledge of the 
model for neutron reflection. 

In all earlier studies in this work, we account for gravity and spatial 
variation in the cold beam fluence. Failure to account for gravity shifts the 
estimate of r* downward by 5.1 s (standard uncertainty is 1.6 s) (Table 4). 
We expect this shift for two reasons. First, gravity changes the distribution 
of the initial potential energy of UCNs and hence the distribution of the 
initial energy of UCNs. Second, for the cases studied here, gravity reduces 
the minimum potential energy on the nominal trapping volume boundary 
from 139 neV to 135.5 neV. Recall, this minimum potential energy is the 
threshold for defining a marginally trapped UCN. 

Failure to account for spatial variation in assumed neutron beam fluence 
shifts the estimate of r* upward by 7.5 s (standard uncertainty is 1.6 s) (Table 
4). As a summary, we list estimates of systematic uncertainties associated 
with the effects studied here in Table 5. 

For the NIST experiment, we expect the walls of the trap to vibrate and 
perturb the energy of UCNs that survive wall collisions and are reflected back 
into the trapping volume. For more discussion of this effect for a simplified 
1-D vibrational model for a magneto-gravitational trap, see Salvat and Wal- 
strom [29]. Since the wall loss probability of an UCN depends on its Ej_ 
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Table 2. For a static potential experiment, dependence of estimated r* on 
lower axial boundary of nominal trapping volume ( z m i n ) and time-step d t in 
symplectic integration algorithm. We account for wall losses when estimating 
t*. Beam divergence effects neglected. Reported uncertainties account for 
uncertainty of A (t) (Eq. 16) due to Monte Carlo sampling variability. 


d t (s) 

Zmin (cm) 

T*(s) 

1.0e-04 

-40 

703.3(1.5) 

1.0e-04 

-35 

700.5(1.5) 

1.0e-04 

-37.5 

701.9(1.1) 

5.0e-05 

-37.5 

699.9(2.4) 

2.5e-05 

-37.5 

699.9(2.8) 


Table 3. For a static potential experiment, dependence of estimated r* on 
neutron scattering modeling. Beam divergence neglected. Reported uncer¬ 
tainties account for uncertainty of A (t) (Eq. 16) due to Monte Carlo sampling 
variability. 


reflection model 

df (s) 

Zmin (cm) 

t*(s) 

diffuse Lambertian 

1.0e-04 

-37.5 

701.9(1.1) 

diffuse uniform 

1.0e-04 

-37.5 

703.2(1.0) 

specular 

1.0e-04 

-37.5 

713.5(1.5) 
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Table 4. For a static potential experiment, dependence of estimated t*. on 
gravity and spatial variation of neutron beam (Figure 2). For all cases, the 
wall loss probability model accounts for incoherent scattering. We assume a 
diffuse Lambertial model for neutron reflections. No beam divergence. The 
time step parameter is df =1.0e-04. The uncertainties account for imprecise 
knowledge of the A (t) (Eq. 16) clue to Monte Carlo sampling variability. 


gravity accounted for 

beam profile accounted for 

r*(s) 

yes 

yes 

701.9(1.1) 

no 

yes 

694.4(1.2) 

yes 

no 

707.0(1.2) 


Table 5. For a static potential experiment, systematic effects and associated 
uncertainties for r* measurement due to wall losses of marginally trapped 
UCNs. Effects that are not statistically significant are indicated with aster¬ 
isks. 


effect 

correction 

uncertainty 

wall loss model 

none 

8 s 

neutron reflections model 

none 

1.3 s 

beam divergence 

-2.9 

2.9 s 

beam profile model 

none 

NA 

upscattering 

none 

<< 1 s 

3 He absorption 

NA 

NA 

mechanical vibrations 

none 

NA 

time step 

none 

1.9 s * 

choice of z mtn 

none 

1 s * 

Total 

-2.9 s 

8.9 s 
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value, energy perturbations would affect the survival probability of above 
threshold UCNs for a static trap experiment. Energy perturbations due to 
mechanical vibrations would also affect how well above threshold UCNs are 
purged (and how well below threshold UCNs are retained) when magnetic 
holds are ramped. Development of a realistic three-dimensional model for 
the effect of mechanical vibrations for the NIST experiment is an important 
and very challenging research topic beyond the scope of this study. 

5. Ramping Studies 

In magnetic trapping experiments, one can purge above threshold UCNs 
from the trap by ramping the trapping potential down and then ramping 
it back up to its original value [4,30]. In a background-free simulation ex¬ 
periment, we reduced the quadrupole held from its maximum value to an 
adjustable fraction of its initial value (Figure 9). We denote this fraction as 
a ramping fraction. For any given ramping schedule, we determined the con¬ 
ditional survival probability for all UCNs (rather than just above threshold 
UCNs) that survive ramping by a direct simulation method for a hxed value 
of r* = 686 s. Based on this conditional survival probability, we simulate 
high-count /3 —decay data under the assumption that the observation rate in 
narrow 1 s width bins is well approximated as (. N(t )) /r n . We ht the Eq. 22 
prediction model to simulated data with A (t) = 0. That is, we ignore wall 
losses. For ramping fractions less than approximately 0.3, the systematic 
error of the estimated lifetime is very small (Figure 10). 

The benefits of purging above threshold UCNs by held ramping come 
with a cost. That is, we also purge a non-negligible fraction of initially below 
threshold UCNs. For instance, for the case where the ramping fraction is 0.3, 
the fraction of surviving above threshold UCNs immediately after ramping 
ends is 0.00072. However, the fraction of originally below threshold UCNs is 
reduced to 0.30(0.02). 

There are other held ramping strategies to purge above threshold UCNs 
than the one considered here. For instance, one could hll the trap with the 
trapping potential reduced to its minimal value and then ramp it after the 
filling stage ends. In a simulation experiment, we compared this alternative 
ramping strategy to ours. We simulated UCNs with the same initial locations 
and initial velocities and tracked them for both strategies. In order to purge 
all but 0.0072 of the above threshold UCNs, the fraction of below threshold 
UCNs retained by the alternative upramping scheme is 0.06(0.01). Thus, 
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for this case, the relative number of below threshold UCNs after ramping is 
less for the alternative ramping strategy compared to the ramping strategy 
implemented in our simulation study. To help explain the result, we note that 
in the alternative ramping schedule, the quadrupole field is maintained at a 
much lower initial value than in the original ramping schedule. Hence, during 
the filling stage, the energy threshold that defines a marginally trapped UCN 
is lower for the alternative ramping scheme compared to the original scheme. 

6. Summary 

We developed a competing risks survival analysis model to account for 
losses due to the joint effect of (3 —decay losses, wall losses of marginal trapped 
neutrons, and an unspecified absorption mechanism. We determined the 
survival probability function associated with the wall loss mechanism by 
a Monte Carlo method based on physical models for the loss probability 
of above threshold UCNs that interact with materials in cylindrical wall 
and endcaps of the trap. In our approach, we track above threshold UCNs 
in the trap with a computer intensive symplectic integration method. For 
any trapping volume, we objectively define above threshold UCNs as those 
with total energy greater than the minimum potential energy on the trap 
boundary. Hence, there is no need to track UCNs below this threshold. 

Based on estimated survival probabilities, we constructed a prediction 
model for observed data acquired in a magnetic trapping neutron lifetime 
experiment at NIST where the trapping potential is static. Based on a fit 
of this model to a subset of the NIST experimental data, we determined the 
mean lifetime of neutrons in the trap to be 700 s with an uncertainty of ap¬ 
proximately 60 s - considerably less than the current best estimate of (880.1 
± 1.1 ) s. The source of this discrepancy is an ongoing research topic; ex¬ 
periments are underway to determine if this discrepancy can be explained by 
neutron capture by 3 He impurities in the trapping volume. In our model, the 
largest source of systematic uncertainty is associated with imperfect knowl¬ 
edge of the wall loss probability model (8 s) (Table 5). If we neglect the 
effect of marginally trapped neutrons, for the static case considered here, the 
estimated lifetime of the neutron is shifted downward by 45 s (Tables 1). 

In a Monte Carlo experiment, we demonstrated that when the trapping 
potential is ramped down and then back again, systematic error due to wall 
losses of marginally trapped UCNs can be suppressed to a very low level. For 
a particular case, Monte Carlo simulation studies showed that this ramping 
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strategy is more efficient than one where neutrons are trapped when the field 
is fully ramped and then increased to its maximum value. From a practi¬ 
cal perspective, our stochastic model and associated Monte Carlo methods 
should be helpful to guide design of field ramping strategies to purge above 
threshold UCNs in magnetic trapping and related trapping experiments. 
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Figure 1: Magnetic field profiles in the NIST magnetic trapping expeirment. 
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Figure 2: Loss probabilities for (a) the Teflon coated endcap; (b) the acrylic coated 
endcap and (c) multilayers at cylindrical wall. 
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Figure 3: Model for neutron fluence in NIST magnetic trapping experiment at trap en¬ 
trance. 
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Figure 4: For a static trap where the trapping volume is defined to be -37.5 cm < z < 
37.5 cm, the ratio of above threshold to below theshold UCNs is 4.22 for energies between 
approximately 139 neV to 246 neV. 
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loss time (s) 



Figure 5: Loss time and wall collision rates' ror simulated above threshold UCNs. A UCN 
is lost when its empirical survival probability falls below 1.0e-09. Tracking halted 5500 s 
after creation time if a UCN has survival probability above 1.0e-09. 
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Figure 6: (a). Monte Carlo estimate of survival probability function of above theshold 

UCNs. (b). Monte Carlo estimate of distortion term A (t) (Eq. 16). 
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Figure 7: Observed and predicted (line) background-corrected data for a subset of data 
from the NIST magnetic trapping experiment. For this subset, the trapping potential is 
static. Prediction model accounts for wall losses. 
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Figure 8: Model Carlo estimates of A (t) depend on assumed beam divergence witin the 
trap. Fitted intercept and slope are 699.73(1.03) and —0.97(0.17) s / degree. 
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Figure 9: In simulation experiment, the quadrupole field is reduced by a fraction that 
varies from 1 to an adjustable minimum. 
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Figure 10: Bias of neutron lifetime when an exponential model is fit to simulated data 
contaminated by above threshold UCNs. For unramped field, bias is —31(2) s. The true 
value of t* is 686 s. 
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